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ABSTRACT 

The last decade has seen applications of Adaptive Mesh Refinement (AMR) methods for a wide range 
of problems from space physics to cosmology. With the advent of these methods, in which space is 
discretized into a mesh of many individual cubic elements, the contemporary analog of the extensively 
studied line radiative transfer (RT) in a semi-infinite slab is that of RT in a cube. In this study we 
provide an approximate solution of the RT equation, as well as analytic expressions for the probability 
distribution functions (pdfs) of the properties of photons emerging from a cube, and compare them with 
the corresponding slab problem. These pdfs can be used to perform fast resonant-line RT in optically 
thick AMR cells where, otherwise, it could take unrealistically long times to transfer even a handful of 
photons. 

Subject headings: line: formation — radiative transfer — resonant 
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1. INTRODUCTION 

The classic problem of resonant radiative transfer (RT) 
in a semi-infinite slab has been extensively studied in liter- 
ature (e.g., Harrington 1973; Neufeld 1990). On the other 
hand, the last decade has seen applications of Adaptive 
Mesh Refinement (AMR) methods to problems as diverse 
as solar physics, supernovae and nucleosynthesis, interstel- 
lar medium physics, star formation, astrophysical jets, cos- 
mology, etc. (for a summary see, e.g., Norman 2004). In 
mesh-based methods the continuous domain of interest is 
discretized into a grid of many individual cubic elements. 
With the advent of AMR codes, the contemporary ana- 
logue - at least in terms of usefulness and applicability 
range - of the extensively studied problem of resonant- 
line RT in a slab is the relatively unexplored problem of 
RT in a cube. 

Understanding resonant RT in optically thick cubes is 
useful in particular because to perform RT in AMR simu- 
lations one must solve numerous cube RT problems, since 
each time a photon enters an AMR cell one has a new 
cube RT problem. Furthermore, as is the case, e.g., for 
Ly-a line RT in cosmological simulations of galaxy forma- 
tion, the high resolution achieved with AMR codes (along 
with the cooling of the gas), leads to very high scatter- 
ing optical depths (see Tasitsiomi 2005). To obtain results 
within realistic times, we need a very fast RT algorithm, 
faster than the standard direct Monte Carlo approach in 
which one follows the detailed scattering of photons in 
each one of the cells. A way to obtain this very fast al- 
gorithm is to study a priori the RT in cubes of various 
physical conditions. Using the results of such a study, we 
can avoid following the detailed photon scattering in each 
cell. Instead, we can immediately take the photon out of 
the cell treating thus the problem on a per cell rather than 
on a per scattering basis, thus accelerating the RT scheme 
considerably. 

In this letter we discuss results on the resonant RT in 
a cube, and in comparison with resonant RT in a slab of 
similar physical conditions. 
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2. DESCRIPTION OF THE PROBLEM AND DEFINITIONS 

In what follows we assume slab and cube configura- 
tions for the scattering medium. The slab is semi-infinite, 
namely it is finite in one spatial dimension and infinite in 
the other two. The scattering medium is uniform in its 
properties and has a central source of center-of-line pho- 
tons that scatter resonantly before escaping. Harrington 
(1973) has solved the slab RT equation in the limit of 
high optical depths and obtained a one parameter solu- 
tion (also, see Neufeld 1990). The parameter is aro, with 
a = Ah , T,/2Ai>r} and Az/j,, Avd the line natural and ther- 
mal Doppler widths, respectively, and To the center-of-line 
optical depth from the center of the scattering medium to 
one of its edges. More specifically, tq is defined so that 
the optical depth at frequency shift Xt — {y — vo)/Avr) 
is T Xf = To0(x/), with 4>{xf) the normalized line profile. 
The discussion that follows applies to optically thick media 
(ar > 10 3 ). 

3. RESULTS 

3.1. Emerging frequency distribution 

3.1.1. Approximate analytic solution for resonant RT in 

a cube 

Following Harrington (1973), one can show that the RT 
equation we need to solve is 

V^(t;<t) + S = -3^(x / )4? (1) 
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where J is the zeroth moment of the intensity I, a is 
defined through dxj/da = (3/2) 1 / 2 ^(a;/), r is defined 
through dr Xf = 4>{xf)dr, j(f) is the emissivity, and f 
is measured from the center of the cube. This equation 
is identical to the equation used previously for a semi- 
infinite slab (Unno 1955; Harrington 1973; Neufeld 1990) 
or a spherically symmetric distribution (e.g., Dijkstra et al. 
2005). The only difference is that those previous cases were 
for one spatial dimension, hence V 2 . consisted of only one 
term, whereas in the case of a cube it has contributions 
from all three dimensions, i.e. V 2 = d 2 J / dr 2 + d 2 J / dr 2 + 
d 2 J/dr 2 , with t x ,t v and t z the components of f along the 
x, y and z-axis, respectively. 2 

2 Thus, for dr = —nods, with s measured along the photon prop- 
agation direction, and kq such that kq4>(xj) is the 'absorption co- 
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Equation (1) is a linear, inhomogeneous, partial differ- 
ential equation. To solve it we use the eigenfunction ex- 
pansion method. Namely, motivated by the method of 
separation of variables (applicable in the case of the cor- 
responding homogeneous equation), we assume that the 
solution can be written as 



J(r;a) = 



CO 

E 

a,/3,7— 1 



X a (T x )Y (T y )Z 7 (r z )G al3l (a) . (2) 



When applying this expansion method in inhomogeneous 
problems, the idea is that X a ,Yjj and Z~ t will be known 
(eigcn)functions and G a pj will be the unknown coefficients 
of the sum (here frequency-dependent) that are to be de- 
termined through the solution process. In reality, we have 
to specify all four functions since we do not have any 
"known" eigenfunctions. We take the "known" (eigen)functions 
of position to be the solutions to the "associated" homo- 
geneous ordinary differential equations. The term "asso- 
ciated" here implies that the best choice for the basis po- 
sition (cigcn)functions should be the solution sets from 
Sturm-Liouvillc problems that closely resemble the prob- 
lem being addressed. This will give a set of functions that 
are orthogonal over the domain defined by the problem. 
Thus, to find the position (eigen)functions, we focus first 
on the solution of the corresponding homogeneous equa- 
tion. We assume that the solution is separable, namely 
that J{f,a) = R(t)G(cf). Substituting in equation (1) 
and after some rearrangement we get 

Vfo(f) 1 d 2 G(a) 2 

R(f) G{&) da 2 ' [ ' 

where we have performed a first separation with separa- 
tion constant —A 2 . As already implied by equation (2), 
we furthermore assume that R(f) = X(t x )Y ~(t v )Z(t z ) in 
Cartesian coordinates. After performing additional sepa- 
rations we end up with equations of the form 

l92X * (4) 



two stream approximation. Now, however, these condi- 
tions refer separately to each one of the three independent 
directions. Thus, we get 

2 dJ 

J(±T ,T y ,T z ;a) = ±2H(±T ,T y ,T z ;a) = TttT/ — r 

6(/>(Xf) OT x 

with H the first moment of / and where, in order to obtain 
the last equality, we used the other Eddington condition, 
J = 3Tr(K), with K the second moment of / with respect 
to n, combined with that V ' T K / '<p(xf) — H. Using equa- 
tion (5) we get for l a the condition (also see, e.g., Neufcld 
1990) 



l a r Q ~ ir(a - 1/2) <^ 1 



3(f)T 



O 



1 



U<M 2 



(6) 



X dr 2 

and similarly for Y(t v ), Z(t z ) with separation constants 
— to 2 and — n 2 , respectively, and with A 2 = n 2 + I 2 + m 2 . 
The general solution for each one of these equations con- 
sists of both sine and cosine terms. Since we only consider 
central point sources, i.e. j(r) = 8{t x )8{t v )5{t z ), each 
one of the functions X,Y, and Z must be separately even. 
Thus, we set X(t x ) = Acos(It x ),Y(t v ) — Bcos(mT y ) and 
Z(t z ) = Ccos(nT z ). 

We calculate the constants A, B, C using boundary con- 
ditions that are generated assuming the Eddington ap- 
proximation (see, e.g., p 322 of Eddington 1926), where 
near isotropy is assumed. Given the optical depths we 
are concern with, the near isotropy assumption should be 
fairly accurate. In fact, we adapt the usual two stream 
approximation to the cube problem. Instead of assuming 
mild anisotropy in the form of two streams (i.e., I = I\(t) 
for < 6 < 7r/2 and I = I 2 for 7r/2 < 6 < n), we take 
into account all 3 spatial directions and treat them equiv- 
alently. That is, we use an eight-stream approximation, 
which nevertheless, leads to the same conditions as the 

efficient', T x ,r y and t z are such that dr x = —Kodx,dr y = —Kody 
and dr z = —Kodz, respectively. In simple terms, t x ,t v and t z are 
(modulo a factor kq that depends on the density of the scattering 
medium) equivalent to spatial coordinates within the cube. 



as well as similar conditions for mp and n 1 . 

Normalizing the cosine solutions in [—To, t ] in each di- 
rection separately, the solution to the homogeneous equa- 
tion RT equation takes the form 



cos(/ Q Ta;) cos(mpT y ) cos(n 7 T z ) 



J(T x ,T y ,T z ;a) = 

a,/3,7— 1 

x G Q/3 »/r 3/2 . (7) 

Now we proceed to calculate the " coefficients" of the sum 
given in equation (2). Substituting this in equation (1), 

and after multiplying by cos(/ ct r :E )cos(m ( 3Tj / )cos(n 7 T2)/TQ^ 2 , 
and integrating over volume we get 

d 2 G n 



-X^^Gafjj 



8a 2 



4ttt, 



(8) 



with 



Qaf3 7 = / / / j(T)cOs(l a T x ) COs(m/3T y ) 

J J —To J —To 

cos(n 7 T z )dT x dTyd,T z . (9) 

For a source of unit strength at the center of the cube 
we have j(r) = 5(T x )5(T y )S(T z ), and thus Qap^ = 1- In 
equation (8) we set 3</>(x/) 2 = \/6S(a) (see, Harrington 
1973). 

Away from a = equation (8) is homogeneous. Impos- 
ing the boundary condition J(r, ±oo) = that reflects the 
fact that we do not expect photons with infinite frequency 
shifts (thus, Gafj-f must go to zero for high - positive or 
negative - a values) we get 

G Q/3 » = De- A «o>l . (10) 

We obtain the value of the constant D exactly as in Har- 
rington (1973). Plugging all these in equation (7) and, 
since we are interested in the overall spectrum of radia- 
tion emerging from one side of the cube, say along the z 
axis, after integrating over t x ,t v and setting t z — To we 
get 



J(t ,<7) 



sm(l a T ) sm(mf3T ) 



a,/3,7— 1 



m f3 T 



cos (n 7 T ) 



y/g g— \J {l a To) 2 + {m«Tn) 2 + {n~,Ta) 2 \o\/Ta 

27r \/(^r ) 2 + (mpTo) 2 + (n 7 r ) 2 



(11) 



At this last step we also substituted A 2 ^ with l 2 a +m 2 ^+n 2 . 
For a comparison of the spectrum emerging from one side 
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of a cube to that emerging from one of the two 'sides' of 
a slab we must multiply our cube solution by a factor of 3 
so that both solutions have the same normalization. This 
is so because, for the same central source, we expect 1/6 
of the photons to emerge from a certain cube side and 1/2 
of the photons to emerge from a certain slab 'side'. 

Each individual series term for both the cube and the 
slab solution (given in Harrington 1973) depends only on 
qt . This becomes obvious when one takes into account 
the definition of a and the approximation for <fi(xf) ~ 
a/nx^ in the wings, as well as condition (6). The slab 
solution is an alternate series and can be written in closed 
form. The cube solution cannot be written in closed form, 
but, using equations (5) and (6) we see that sin(^ a ro) ~ 
(-1)"- 1 and cos(n 7 r ) = 2n 7 /30(x / )(-l)''- 1 . Thus, the 
sin(7 a To) sin(m ( 3To) cos(n 7 ro) term in equation (11) can be 
written as — (— l) a +P+T, indicating that the cube series 
may also be alternating. Writing this three-sum series as 
one sum, i.e., YmLi c i-> we nn d that indeed the cube series 
is alternating as well. Some of the first finite sums of the 
alternating series for a cxtq = 2 x 10 4 cube and slab are 
shown at the top two panels, and the left bottom panel of 
Figure 1, respectively. In the case of the cube the infinite 
number of terms result (solid histograms) is obtained by 
the Monte Carlo method described in detail in Tasitsiomi 
(2005), whereas for the slab we use the closed form slab 
analytical solution of Harrington (1973). 

The cube series solution will be of some practical value, 
only if a few first terms contribute significantly to the 
sum. The exponential decay in \a\ indicates that the terms 
should die off for 'high' a, (3 and 7 values-with the ex- 
act values where this happens dependent on the frequency 
(or a) one calculates the spectrum at. The logarithm of 
the absolute ratio of the series terms, Cj, in units of the 
first term, c\ , for 3 different values of the frequency shift 
is shown in the bottom right panel of Figure 1. As be- 
fore, ar = 2 x 10 4 , but we find that these 'convergence 
curves' remain identical for all cube aro at the extremely 
optically thick regime. The dashed-dotted line is the 'con- 
vergence curve' for a slab for frequency shift equal to the 
shift where the emergent spectrum is known to have a 
maximum (<~ 0.9(ar ) 1 / 3 ; Harrington 1973). Away from 
this maximum, the slab convergence curves look identical 
to those of the cube (dotted and dashed lines). Clearly, 
the rate of convergence depends on the frequency the spec- 
trum is calculated at. Given that the series is alternating, 
with the absolute value of the terms decreasing monoton- 
ically (for most a values), the absolute error we make by 
truncating the series at the ith term is less or equal to 
the i+1 term. From the plot, this means that if we only 
keep the first term, the actual infinite series limit S will be 
~ (lilO- 4 )^! for x f = 39 and - (l±0.3)Si for x = x maxi 
with Si the ith partial sum. We can achieve a 3% accu- 
racy at the peak if we go to i = 4 (i.e., \S — S^/S < 0.03), 
whereas for better accuracy we have to exceed i = 30. 
Convergence becomes extremely slow at values close to 
the resonance, e.g., x = 2. Furthermore, we find that the 
up to i = 30 terms for x = 2 do not decrease monoton- 
ically in the case of the cube - they do decrease for the 
slab but extremely slowly. 

The usefulness of the slab series solution derived here 
depends on the application in mind. For example, when 



studying Ly-a emitters at very high redshifts, it is ex- 
pected that the blue frequencies will be absorbed anyway 
because of hydrogen intervening between the emitter and 
the observer, and the red wavelengths near the resonance 
will also be absorbed by the red damping wing. In this 
case, the poor convergence at frequencies near resonance 
(and, more generally, over the |a;| < x max range) may not 
be a problem. Clearly, one must decide on the usefulness 
of the series solution based on the specifics of the applica- 
tion. 
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Fig. 1. — Top left panel: Cube emergent spectrum obtained via 
Monte Carlo RT simulation (solid histograms, = 00), and approxi- 
mate analytical spectra using the one or two first terms of the sum 
given by equation (11) (dotted and dashed line, respectively). The 
n = 1 spectrum is multiplied by 0.5 to fit within the same y-axis 
range as the other spectra. Top right panel: Same as in left panel 
but for n = 30 terms (dashed line). Also shown is the spectrum 
obtained using the slab solution for a smaller, effective aro (dotted 
line). Bottom left panel: Slab emergent spectra using for the slab 
infinite number of terms (solid line— obtained using the closed form 
solution), one term (dotted line), two terms (dashed line), and the 
first thirty terms (dot-dashed line). Bottom right panel: Logarithm 
of the absolute value of the ith term of the sum representing the 
spectrum for a cube/slab, in units of the first term ci. Results are 
shown for the spectrum estimated at different frequency shifts, x. 
For most x values we show results only for the cube, since cube and 
slab terms exhibit the same convergence behavior, except near the 
frequencies where they attain their maximum intensity (x m ax)- All 
results are for qtq = 2 X 10 4 . See text for details. 



3.1.2. Approximate cube spectrum in closed form 

Based on the similarities between the cube and the slab 
spectra, one might think that the spectrum emerging from 
one side of a cube may be well described by the slab closed 
form solution, but for a different, smaller aro than the 
actual ar of the cube. Because, for example, when ob- 
serving the spectrum emerging from a cube along the z- 
direction we lose all photons that in the case of a slab 
would wander, scatter many times along the infinite di- 
mensions until they finally escape from the z-plane. In 
the case of the cube these photons do not contribute to 
the spectrum we observe from the z-direction because; they 
have already exited the cube along directions other than z. 
This argument has been already used in Tasitsiomi (2005), 



4 



TASITSIOMI 



where the slab solution was used to describe the spectrum 
for a cube by plugging in the slab solution an ar equal 
to 2/3 the actual cube ar . This value was motivated by 
comparing the mean number of scatterings, N sc , for pho- 
ton escape in cubes and slabs. Since N sc scales linearly 
with t for extremely optically thick media (Adams 1972), 
by comparing N sc for cubes and slabs one can guess a cor- 
rect effective r (and thus aro). For the purposes of this 
previous study, this effective thickness gave good agree- 
ment with simulated spectra for a wide range of physical 
conditions. 

Attempting a more detailed treatment, we fit cube spec- 
tra with the slab solution for the fraction / of aro that 
must be used in this solution to get the best fit. We find 
that / varies in the 0.66 — 0.77 range for aro values in 
the 2 x 10 3 — 10 8 range. We find no trend of / with ar . 
Within errors, one can use a constant fraction in the above 
/-range regardless of aro since we find that we cannot dis- 
tinguish between the fits obtained when 2/3 or the exact / 
value are used. An example of a fit of the cube spectrum 
using the slab solution and / = 2/3 is shown at the top, 
right panel of Figure 1 (the best fit /-value for this exam- 
ple is 0.72). Furthermore, in terms of the spectrum shape 
(e.g., maximum location and intensity), the above range 
of / values lead to spectra indistinguishable by currently 
existing instruments. 

3.2. Distribution of exiting direction and point 

Referring to /i, the cosine of the angle with which the 
photon is exiting, measured with respect to the normal 
to the exiting surface, we show its cumulative probability 
distribution function (cpdf) at the top panel of Figure 2. 
We find that this cpdf is the same for a slab or a cube, 
and clearly deviates from isotropy (dashed line) . We verify 
the findings of other studies that in optically thick media 
photons tend to exit in directions perpendicular to the 
exiting surface (see, e.g., Phillips & Meszaros 1986; Ahn 
et al. 2002). 

In cases of very optically thick media, the emerging radi- 
ation directionality approaches the Thomson scattered ra- 
diation emergent from a Thomson-thick electron medium. 
Phillips & Meszaros (1986) found that for a Thomson- 
thick medium I(p)/I{0) = (1 + 2ju)/3. Since the number 
of photons emerging at \i is cx I{pL)[id\i, we get 

^(<M)= £ + o 2M w" = Y( 3 + 4 ^- < 12 > 
J (1 + 2/i)/xd/x 7 

This dependence is shown in Figure 2 with the solid line. 
It is an excellent description of the directionality of the 
spectrum and hence equation (12) can be used to deter- 
mine the photon exiting direction. It has been suggested 
by some authors (Ahn ct al. 2002) that the fact that in 
optically thick media RT occurs mostly via wing photons 
with the latter being described by a dipole phase function 
(Stenflo 1980), and the fact that Thomson scattering is 
also described by a dipole scattering phase function ex- 
plain why the resulting [i cpdfs are similar. However, we 
find the same cpdf if we use either an isotropic or a dipole 
phase distribution. For such optical thicknesses the de- 
tails of the exact phase function do not matter, at least 
not with respect to the exiting angle cpdf. Regardless of 
scattering phase function, in such media photons tend to 




Fig. 2. — Top panel: Cumulative probability distribution func- 
tion of the cosine of the angle the photon exiting direction forms 
with the normal to the exiting surface (fi). Results arc shown for 
a semi-infinite slab or a cube (long- dashed and short- dashed lines, 
respectively), as well as for an isotropic or dipole scattering phase 
function (long-dashed and dot-dashed line, respectively). With the 
solid line we show the directionality of radiation emerging from a 
Thomson-thick electron medium. All these distributions arc indis- 
tinguishable. Also shown is the case where photons exit the slab (or 
cube) isotropically (dashed line). Bottom panel: Cumulative proba- 
bility distribution function of exit points of radiation emerging from 
a slab (dashed line) and a cube (dotted line) of the same aro- The 
photons are assumed to emerge along the z direction, hence the dis- 
tribution refers to either the x or y coordinate, in units of the size 
of the cube (or the finite dimension of the slab). 



escape along the normal to the exiting surface where the 
opacity is smaller. The results shown are for ar = 2 x 10 3 , 
but similar distributions are obtained for thicker media. 

In the case of the slab the azimuthal angle 4> is dis- 
tributed uniformly in [0, 2ir]. In the case of the cube there 
are small deviations from uniformity. This is expected 
since the previous distribution for fi is found to be valid 
for all sides of the cube. In the simplest case where we 
observe along the z-axis (in which case the spherical coor- 
dinate <f> angle is the actual azimuthal angle we refer to), all 
direction cosines - cos 9 and sin 9 cos <f> and sin 9 sin <j) - fol- 
low the distribution given in equation (12), thus <j> cannot 
be exactly uniformly distributed. However, the deviations 
from uniformity are small. 

The distribution of exiting points for both a cube and a 
slab are shown at the bottom panel of Figure 2. For this 
figure we assume that we observe photons emerging along 
the z direction and we record the x and y coordinates of 
their exit points (in units of the size of the cube). In the 
case of RT in a cube the distribution is pretty close to 
uniform. In the case of the slab, despite it being semi- 
infinite, for any practical purposes one can assume that 
all photons exit at most within |a;|(|y|) < 5 (not shown in 
figure). 
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